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Abstract.  Three  dimensional  (3D)  FLASH  Laser  Radar  (LADAR)  sensors  are  unique  due  to 
the  ability  to  rapidly  acquire  a  series  of  two  dimensional  remote  scene  data  (i.e.  range  images). 
Principal  causes  of  3D  FLASH  LADAR  range  estimation  error  include  spatial  blur,  detector 
blurring,  noise,  timing  jitter,  and  inter-sample  targets.  Unlike  previous  research,  this  paper  ac¬ 
counts  for  pixel  coupling  by  defining  the  range  image  mathematical  model  as  a  2D  convolution 
between  the  system  spatial  impulse  response  and  the  object  (target  or  remote  scene)  at  a  partic¬ 
ular  point  in  time.  Using  this  model,  improved  range  estimation  is  possible  by  object  restoration 
from  the  data  observations.  Object  estimation  is  performed  by  deriving  a  blind  deconvolution 
Generalized  Expectation  Maximization  (GEM)  algorithm  with  the  range  determined  from  the 
estimated  object  by  a  normalized  correlation  method.  Theoretical  derivations  and  simulation 
results  are  verified  with  experimental  data  of  a  bar  target  taken  from  a  3D  FLASH  LADAR 
system  in  a  laboratory  environment.  Simulation  examples  show  that  the  GEM  improves  range 
estimation  over  the  unprocessed  data  and  a  Wiener  filter  method  by  75%  and  26%  respectively. 
In  the  laboratory  experiment,  the  GEM  improves  range  estimation  by  34%  and  18%  over  the 
unprocessed  data  and  Wiener  filter  method  respectively. 

Keywords:  3D  FLASH  LADAR,  laser  radar,  range  estimation,  blind  deconvolution,  object 
restoration,  waveform  processing,  generalized  expectation  maximization 

1  INTRODUCTION 

A  three-dimensional  (3D)  FLASH  laser  radar  (LADAR)  is  a  pulsed  radar  system  this  is  both 
an  imaging  and  ranging  sensor.  Referring  to  Fig.  1,  a  3D  FLASH  LADAR  produces  a  time 
sequence  of  two-dimensional  (2D)  images  due  to  a  fast  range  gate  resulting  in  a  3D  data  cube 
of  spatial  and  range  scene  data  with  excellent  range  resolution  [1],  [2],  FLASH  technology 
principally  differs  from  scanning  LADAR  by  being  able  to  form  a  3D  representation  of  a  re¬ 
mote  scene  in  one  laser  pulse  rather  than  rastering  a  3D  scene  together  using  many  pulses. 
This  capability  results  in  faster  scene  collection  times  with  lighter  weight,  lower  power,  and 
reduced  mechanical  complexity  as  compared  to  the  scanning  systems.  Practical  applications  of 
3D  FLASH  LADAR  include  intelligence,  surveillance,  and  reconnaissance  (ISR),  rendezvous 
and  capture  (air  and  space),  collision  avoidance  (COLA)  for  manned  and  unmanned  air/ground 
vehicles,  and  weapon  targeting. 

Typically,  a  3D  FLASH  LADAR  operates  in  one  of  two  modes.  The  first  mode  is  called 
“HIT  mode”  where  each  pixel  element  (pixel)  is  independently  triggered  when  its  intensity 
reaches  a  preset  threshold.  This  mode  is  advantageous  when  searching  for  a  target  where  the 
range  is  not  already  known.  However,  truncated  waveforms  can  occur  leading  to  range  estima¬ 
tion  errors.  The  second  mode  is  called  “SULAR  mode”  where  the  pixels  are  triggered  to  start 
recording  data  together  based  on  a  preset  range.  Benefits  of  this  mode  include  being  able  to  suc¬ 
cessively  capture  fine  details  of  the  target  and  background.  Drawbacks  are  that  the  target  range 
must  be  known  a  priori  and  waveforms  are  truncated  for  targets  near  the  end  of  the  collect. 
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Fig.  1.  (a)  3D  view  of  LADAR  system  model  in  Cartesian  coordinates  with  each  data  cube  having 
dimensions  of  pixel  x  pixel  x  time  sample.  The  variable  d(tk)  corresponds  to  the  kth  receiver  detected 
range  image  with  k  £  [1,  N],  (b)  Another  view  of  the  3D  FLASH  LADAR  operation.  Each  range  image’s 
full  field  of  view  (FOV)  is  128  x  128  pixels  with  a  range  gate  near  2  nanoseconds  corresponding  to  the 
3D  FLASH  LADAR  system  used  for  experimental  collects. 
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A  method  to  model  the  3D  FLASH  LADAR  data  operating  in  SULAR  mode  is  that  the 
2D  range  images  are  formed  via  a  convolution  between  the  object  at  a  particular  time  and 
the  spatial  impulse  response.  In  Fig.  1(a),  a  range  image  d(tk)  is  one  of  the  2D  slices  of  the 
data  cube.  Considering  the  laser  illuminating  a  target,  one  collect  from  a  3D  FLASH  LADAR 
sensor  results  in  a  data  cube  consisting  of  a  series  of  range  images  (N  from  Fig.  1)  representing 
detected  photons. 

Attempts  at  3D  FLASH  LADAR  range  estimation  of  a  remote  scene  can  result  in  errors  due 
to  several  factors  including  the  optical  spatial  impulse  response,  detector  blurring,  photon  noise, 
timing  jitter,  and  readout  noise.  These  factors  either  cause  the  scenes  intensity  to  spread  across 
pixels  or  add  unwanted  and  disruptive  noise  effects.  The  intensity  spreading  and  noise  corrupts 
the  correct  pixel  intensities  by  mixing  intensities  with  neighboring  pixels  thereby  providing 
false  intensity  values  and  therefore  incorrect  photon  counts  to  the  range  estimator.  Without  blur 
and  noise  compensation,  the  range  estimates  would  be  inaccurate  to  a  degree  depending  on  the 
blur  and  noise  severity. 

The  motivation  behind  this  paper  is  to  provide  a  means  of  improving  range  estimation  by 
object  recovery  (i.e.  spatially  deblurring  data)  from  3D  FLASH  LADAR  observations.  Refer¬ 
ring  to  Fig.  1(a),  the  idea  is  to  process  the  data  in  the  spatial  dimensions  ( x ,  y)  while  improving 
ranging  performance  in  the  time  dimension  ( z ). 

The  theoretical  development  of  the  range  estimator  algorithm  is  covered  first  and  then  ver¬ 
ified  using  simulation  and  experimental  results.  The  algorithm  is  a  variation  on  the  Expecta¬ 
tion  Maximization  (EM)  algorithm  called  Generalized  Expectation  Maximization  (GEM)  and 
is  desirable  due  to  its  iterative  likelihood  maximization,  convergence  properties,  and  ability  to 
decouple  terms  [3].  The  GEM  algorithm  is  powerful  in  that  it  can  perform  blind  deconvolu¬ 
tion  in  situations  with  severe  defocus  or  other  aberrations  including  atmospheric  turbulence. 
To  account  for  different  scenarios,  two  versions  of  the  GEM  algorithm  are  derived  that  either 
recover  the  pulse-shape  or  the  object.  The  primary  difference  between  the  two  involves  data  re¬ 
quired  and  accuracy.  Pulse-shape  estimation  requires  less  data,  but  is  less  accurate  than  object 
estimation.  Additional  details  of  the  differences  are  presented  in  Sections  2.3  and  2.4. 

In  addition  to  the  GEM  algorithms,  a  Wiener  filter  method  is  used  to  attempt  range  estima¬ 
tion  improvement  via  object  recovery  from  3D  FLASH  LADAR  observations  [4],  [5].  Requir¬ 
ing  spatial  impulse  response  knowledge  a  priori,  this  method  can  only  perform  deconvolution 
unlike  the  blind  deconvolution  ability  of  the  GEM.  The  purpose  for  adding  this  other  method  is 
to  show  that  the  GEM  outperforms  a  competing  algorithm  that  already  knows  part  of  the  answer 
(spatial  impulse  response). 

2  THEORETICAL 

This  section  discusses  incoherent  imaging  and  the  application  to  3D  FLASH  LADAR,  presents 
the  data  model,  and  develops  the  GEM  algorithms  that  perform  object  recovery  leading  to  im¬ 
proved  range  estimation.  Even  though  the  laser  light  is  partially  coherent,  the  argument  is  made 
that  the  detected  light  is  able  to  be  modeled  as  fully  incoherent.  Consequently,  this  allows 
for  the  returns  to  be  a  result  of  a  linear,  spatially  invariant  (LSI)  system  involving  an  intensity 
convolution  (instead  of  amplitude  convolution)  between  the  intensity  point  spread  function  and 
the  remote  scene.  Linearity  is  a  consequence  of  electromagnetic  wave  propagation  theory,  and 
spatial  invariance  results  from  remaining  with  the  isoplanatic  angle  [6].  Utilizing  this  LSI  con¬ 
volution  model,  two  GEM  blind  deconvolution  algorithms  are  developed  that  enable  improved 
range  estimation. 

2.1  Incoherent  imaging 

Both  the  negative  binomial  and  Poisson  distributions  can  be  used  to  capture  the  non-negative, 
discrete  nature  of  the  laser  light.  The  negative  binomial  distribution  would  be  the  most  optimal 
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in  describing  the  illuminating  partially  coherent  laser  light,  but  blind  deconvolution  methods  are 
cumbersome  [7].  Whereas,  blind  deconvolution  methods  with  the  Poisson  distribution  (incoher¬ 
ent  imaging)  are  more  tractable  and,  thus,  utilized  in  this  research.  Even  if  the  speckle  is  severe, 
the  benefit  of  modeling  the  speckle  does  not  outweigh  the  cost  of  implementing  a  partially  co¬ 
herent  blind  deconvolution  model  for  the  3D  FLASH  LADAR  system.  Previous  research  using 
the  incoherent  data  model  for  a  3D  FLASH  LADAR  has  also  experienced  success  [8],  [9]. 

To  gain  more  insight  into  this  assumption,  a  simple  approach  is  to  estimate  the  amount  of 
coherence  contained  within  the  3D  FLASH  LADAR  data  by  estimating  the  speckle  parameter 
of  the  negative  binomial  distribution  directly  from  the  data  [7].  Capturing  both  temporal  and 
spatial  coherence,  if  the  speckle  parameter  estimate  is  high  enough,  the  negative  binomial  dis¬ 
tribution  will  look  Poisson-like  allowing  the  data  observations  to  be  modeled  as  arising  from  an 
intensity  convolution  (incoherent  imaging).  Including  speckle  and  photon  noise  effects,  the  neg¬ 
ative  binomial  probability  mass  function  (PMF)  describes  the  photon  distribution  of  a  partially 
coherent  imaging  system  for  a  single  pixel  or  [7] 


P{K) 


T(K  +  M) 
Y(K  +  l)r(7W) 


(1) 


where  M  is  the  speckle  parameter  and  K  is  the  pixel’s  average  photon  count.  Changing  the 
distribution  for  a  3D  FLASH  LADAR,  the  illuminating  laser  light  statistics  for  a  particular 
volume  element  (voxel)  (x.  y,  k  remain  constant)  across  many  data  cubes  is 


P(Djk  ( x,y )  =  djk  (x,  y)  Mj  <E  (1,2, ...,  J))  = 


tt  T  ( djk  ( x ,  y)  +  M) 
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where  j  represents  the  data  cubes,  k  is  the  range  image  (i.e.  time  variable)  within  a  data  cube, 
(x,  y)  are  the  coordinates  in  the  image  plane,  and  djk  (x,  y )  is  the  data  observation.  The  voxels 
are  assumed  statistically  independent  from  each  other  because  of  the  discrete  nature  of  pho¬ 
tons  and  the  detected  photons  do  not  affect  future  detected  photons.  The  maximum  likelihood 
solution  for  the  average  voxel  intensity  is  determined  by 


K  =  7 (x,y). 
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Taking  the  natural  log  of  Equation  (2)  yields 
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where  graphical  methods  are  employed  to  find  the  speckle  parameter  that  maximizes  this  log- 
likelihood.  Using  the  same  experimental  data  as  in  the  range  estimation  efforts,  a  collection 
of  voxels  with  the  strongest  laser  light  are  chosen  to  estimate  the  speckle  parameter.  Figure  2 
shows  the  similarities  between  the  negative  binomial  and  Poisson  distribution  using  an  average 
of  the  estimated  speckle  parameter. 

Even  without  considering  speckle  parameter  estimation  results,  the  argument  can  be  made 
for  incoherent  imaging  due  to  the  Poisson  distribution’s  ability  to  model  the  non-negativity  and 
discrete  nature  of  light  [8],  [9].  This  argument  is  solidified  by  the  speckle  parameter  estimation 
results  indicating  that  the  speckle  noise  appears  low  enough  for  the  incoherent  imaging  model 
to  be  used  with  confidence. 
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Fig.  2.  This  plot  shows  the  negative  binomial  (NB)  using  an  estimated  average  speckle  parameter  ver¬ 
sus  the  Poisson  distribution  with  the  same  mean.  While  not  identical,  the  negative  binomial  distribution 
compares  well  enough  to  the  Poisson  distribution  to  assume  incoherent  imaging. 

2.2  Data  model 

Considering  a  3D  FLASH  LADAR  sensor  with  statistically  independent  samples,  the  PMF  of 
the  observed  photons,  djk  (x,  y),  incorporating  all  cubes  (j  £  [1,  J]),  range  samples  (k  £ 
[1  ,K]),  and  detector  pixels  (x  £  [l,2f],y£  [l,F])is 

P  [Djk  (x,  y)  =  djk  (x,  y) ;  Vj,  k,  x,  y]  = 

TT  [■ ijk  (x,  y)  +  B  (x,  y)]d^x'v)  exp  [ijk  (x,  y)  +  B  (x,  y)]} 

iiw  ^*  (**»)' 

where  the  data  observations  are  defined  by 

djk  (x,  y)  =  ijk  (x,  y)  +  B  (x,  y)  +  njk  (x,  y)  (6) 

with  ijk  (x,y)  as  the  blurry,  non-noisy  data,  B  (x,y)  as  the  pixel  bias,  and  rijk  (x,y)  as  a 
general  noise  term  accounting  for  all  noise  sources  (i.e.  photon  noise,  thermal  noise,  readout 
noise,  etc.).  The  (x,  y)  and  k  variables  correspond  to  a  pixel  in  the  detector  array  and  to  the 
returned  signal  time  of  arrival  respectively.  The  time  of  arrival  is  computed  based  on  the  time 
from  laser  pulse  transmission  to  photon  detection.  The  blurry,  non-noisy  data  is 

M  N 

ijk  (x,  y)  =  °k  (m> n )  ho  {x  -  m,y  -  n)  (7) 

m= 1 n= 1 

where  the  object  ok  (m,  n )  is  defined  at  the  object  plane  with  coordinates  (m  £  [1,  M\ ,  n  £ 
[1,  TV])  and  changes  within  a  single  cube,  but  is  considered  constant  across  all  cubes.  This 
assumption  requires  the  ability  to  perform  cube  registration  due  to  the  possibility  of  moving 
targets,  moving  sensor  platform,  or  inter-cube  timing  errors.  Incorporating  contributions  from 
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Fig.  3.  3D  FLASH  LADAR  data  block  diagram.  For  an  arbitrary  cube  ( j  notation  ignored),  the  observed 
data  is  an  addition  0  of  the  pixel  bias,  the  noise,  and  the  blurry,  non-noisy  image  (convolution,  depicted 
by  0,  between  object  and  PSF). 


light  propagation,  optical  abberations,  and  atmospheric  blurring,  the  intensity  point  spread  func¬ 
tion  (PSF)  hj  (x,  y)  is  constant  within  a  single  cube  while  changing  across  cubes.  In  this  re¬ 
search,  the  PSF  is  considered  constant  within  a  single  cube  since  collection  times  spans  under 
forty  nanoseconds  and  any  time-dependent  effects  would  be  essentially  frozen.  In  addition,  the 
pixel  bias  B  ( x ,  y)  is  constant  between  cubes  as  well  as  within  a  single  cube  due  to  the  pixel’s 
unchanging  physical  material  and  response  to  incident  light. 

Every  pixel  in  the  detector  array  records  a  time-delayed  and  attenuated  version  of  the  trans¬ 
mitted  pulse.  The  physical  outgoing  pulse  shape  of  a  3D  FLASH  LADAR  is  either  Gaussian, 
negative  parabolic,  or  some  hybrid  of  the  two.  The  object  can  be  defined  by  an  amplitude  term 
and  a  pulse  shape  or 

ok(m,n)  =  A(m,n)pk(m,n) .  (8) 

Assuming  a  Gaussian  transmitted  pulse,  the  object  is 
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where  A  (to,  n)  is  the  object  amplitude,  aw  is  the  waveform  standard  deviation,  tk  is  the  time 
variable,  c  is  the  speed  of  light,  and  R  (to,  n)  is  the  range  to  the  target.  For  military  targeting 
or  navigation,  range  to  target  (located  in  the  object  term)  is  the  desired  unknown  variable.  For  a 
negative  parabolic  waveform  shape,  the  object  is  defined  by 
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where  2 pw  is  the  pulse  width  and  rect  is  the  rectangle  function  defined  by 


[0,  if  |x|  >1/2 

rect  (x)  =<  1/2,  if|ir|  =  l/2  (11) 

U  if  M  <  1/2- 

In  attempting  to  perform  range  estimation,  a  range  term  is  not  explicitly  in  the  model,  but  it  is 
buried  within  the  object,  ok  (to,  n),  term  given  by  Equations  (9)  and  (10).  If  the  object  were 
known,  the  target  range  could  be  then  extracted  from  the  object  by  peak  detection  methods. 
Therefore,  the  goal  in  this  paper  is  to  accurately  estimate  the  object  and  recover  the  range 
by  using  a  modified  version  of  peak  detection  that  permits  sub-sample  ranging.  Sub-sample 
ranging  is  done  to  account  for  targets  that  are  between  adjacent  recorded  samples. 

The  unknown  parameters  in  this  estimation  scenario  are  the  object  (target  amplitude  and 
target  range),  PSF,  and  pixel  bias.  The  variable  of  interest  in  this  paper  is  the  range  term  resid¬ 
ing  in  Equation  (9)  or  (10).  Direct  estimation  of  the  range  term  is  problematic  because  of  its 
location  either  in  an  exponential  or  in  a  squared  term.  Therefore,  the  approach  to  range  estima¬ 
tion  is  to  retrieve  the  range  from  the  estimated  pulse-shape  or  object.  This  methodology  relies 
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on  the  knowledge  that  the  target  produces  the  waveform  peak  in  the  detected  returns.  Con¬ 
cerning  the  PSF,  blind  deconvolution  techniques  must  be  employed  since  the  PSF  is  unknown. 
Blind  deconvolution  has  a  rich  heritage  in  astronomical  imaging  providing  a  bevy  of  literature 
attempting  blind  deconvolution.  Although,  blind  deconvolution  in  astronomical  cases  consists 
of  trying  to  recover  one  object  and  one  PSF  (or  many  PSFs  if  using  multiple  frames).  In  trying 
to  recover  the  target  range  from  one  3D  FLASH  LADAR  data  collect,  this  problem  consists  of 
many  objects  with  one  PSF.  There  are  many  objects  due  to  the  transmitted  waveform  causing 
each  range  slice  to  contain  different  intensities  corresponding  to  where  the  waveform  is  incident 
on  the  object.  Therefore,  these  incident  points  become  distinct  objects  in  the  blind  deconvolu¬ 
tion  framework.  If  multiple  cubes  are  necessary,  the  atmosphere  is  changing  with  each  cube 
resulting  in  multiple  PSFs  that  must  be  estimated  resulting  in  a  “multi-frame”  or  “multi-cube” 
scenario.  If  no  atmospheric  turbulence  exists  or  is  non-volatile,  the  PSF  is  consistent  throughout 
the  cubes  and  the  j  subscript  can  be  dropped. 

Traditional  linear  maximum  likelihood  efforts  do  not  suffice  to  estimate  target  range  given 
the  unknowns  in  the  statistical  model  depicted  by  Equation  (5).  More  powerful  object  estima¬ 
tion  methods  like  the  EM  algorithm  must  be  employed  due  to  the  coupled  unknowns  which 
will  be  covered  in  the  next  section.  However,  a  closed  form  solution  for  the  EM  algorithm’s 
maximization  step  is  intractable.  Consequently,  the  GEM  algorithm  goal  is  to  modify  the  EM 
structure  such  that  the  likelihood  is  incrementally  increased  rather  than  globally  maximized  as 
in  the  EM  algorithm.  This  incremental  increase  in  the  likelihood  simplifies  the  maximization 
step  allowing  unknown,  non-random  parameter  estimation.  In  summary,  the  GEM  algorithm  is 
written  as 

Q  >  Q  (V”);  w)  (12) 

where  is  the  vector  of  unknown  variables,  v  is  the  iteration,  and  Q  is  the  expected  value  of 
the  complete  data  log-likelihood  or 

Q  (*;  *(,;))  =  £*(.)  {In Lcd  (<%)}  (13) 

with  Lcd  as  the  complete  data  likelihood  and  y  as  the  incomplete  data.  Complete  data  can  be 
viewed  as  the  unobserved  variables  (fabricated  or  not)  used  to  simplify  the  problem.  Incomplete 
data  is  usually  the  observed  data.  If  Equation  (12)  holds,  it  has  been  shown  that  the  likelihood 
is  increased  with  every  iteration  or  [3] 

L  (Vfe+1))  >  L  (14) 

and,  if  bounded,  the  GEM  sequence  converges  to  a  local  maximum  due  to  the  monotonicity  of 
the  algorithm. 

With  the  target  range  extracted  from  the  estimated  object,  object  recovery  is  accomplished 
using  two  approaches  concerning  the  pulse-shape  and  object  variables  from  Equation  (9).  The 
pulse-shape  estimation  is  very  powerful  in  that  the  estimator  only  needs  one  data  cube  (one- 
shot,  one-kill).  However,  if  the  best  accuracy  is  required  and  the  3D  data  cubes  are  properly 
registered,  the  multi-cube  object  estimation  provides  lower  error. 

2.3  Range  estimation  using  pulse-shape  recovery  via  the  GEM  algorithm 

Considering  pulse-shape  recovery  with  only  one  cube  required  for  processing  (ignoring  j,  since 
j  =  1),  the  problem  is  reformed  by  calling  the  original  data,  dk(x,  y),  the  incomplete  data  and 
specifying 

M  N 

dk  (x,y)  =  EE  dk{x,y\m,n)  +  qk{x,y)  (15) 

m—1  n=  1 
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where  two  new  variables,  dk{ x,  y\m,  n )  and  q k  (x,  y),  are  called  complete  data.  This  formula¬ 
tion  provides  two  sets  of  complete  data  that  account  for  the  photon  noise/image  formation  and 
pixel  bias  respectively.  The  formation  of  the  complete  data  highlights  the  powerful  nature  of 
the  EM  algorithm.  In  this  application,  complete  data  can  also  be  called  unobserved  data  and 
carries  no  explicit  physical  meaning.  It  is  used  to  directly  benefit  the  EM  algorithm.  Consistent 
with  [10],  careful  definition  of  the  complete  data  allows  decoupling  of  unknown  variables  while 
preserving  physical  meaning  in  the  expected  value  of  the  incomplete  data. 

The  expected  value  of  the  complete  data  sets  is  given  by 


E 


dk  ( x,y\m,n ) 


A  (to,  n)  pk  (m,  n)  h  (x  —  m,y  —  n) . 


(16) 


and 

E  [qk  {x,  y)]  =  B  (x,  y)  (17) 

where  B  (x,  y)  is  the  constant  pixel  bias.  The  expected  value  of  the  incomplete  data  is  thus 

E  [dk  (x,  y)]  =  ik  (x,  y)  +  B  (x,  y)  (18) 

which  is  consistent  with  the  data  observations  depicted  in  Fig.  3.  Adding  the  pixel  bias  to  the 
model  covers  non-modeled  noise  effects  and  pixel-to-pixel  impulse  response  variations.  The 
pixel  bias  is  assumed  to  be  governed  by  the  Poisson  distribution  due  to  the  discrete  random 
nature  of  dark  current  and  electron  noise.  Physically,  the  pixel  bias  is  added  to  the  photons 
incident  upon  the  detector  and  is  part  of  the  detected  photon  counts.  The  PMF  for  the  photon 
noise  is 

P  (dk  (x,  y\m,  n)j  = 

[A  (to,  n)  Pk  (to,  n)  h  (x  —  m,y  —  n)]dk{x'y)  eAMm,n)pk(.rn,n)h(x-m,y-n)] 

dk(x,y)\  (19) 


while  the  PMF  for  the  pixel  bias  is 


P(Qk  (x,y)) 


B(x,  yyk(x’y) 
dk  (x,y)\ 


(20) 


Assuming  statistical  independence  between  all  the  pixels  and  between  the  photon  noise  and 
pixel  bias  noise,  the  complete  data  log-likelihood  function  considering  all  random  variables  is 


Lcd  (pk,A,h,B)  =  In 


p(dk(x,y\m,n)^P(qk(x,y)) 


(21) 


k,x,y,m,n 

or  (NOTE:  summations  wrap  around  unless  otherwise  stated) 

Lcd  (pk,A,h,B)  =  ^  dk  (x,  y\m,  n)  In  [A  (m,  n)pk  (to,  n)  h  {x  -  to,  y  -  n)] 

k,x,y,m^n 

-  [A  (to,  n)  pk  (to,  n)h(x  —  m,y  —  n)]  +  qk  {x,  y)  In  [B  (x,  y)]  -  B  (x,  y) .  (22) 
Referring  to  Equation  (13),  the  Q  function  then  becomes 

Q  ( Pk ,  A,  h,B)  =  E  [Lcd  ( Pk,A ,  h,  B)  \dk  (x,  y)  ,p°kld,  Aold ,  hold,  B°ld]  (23) 
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where  the  estimates  for  the  amplitude,  pulse-shape,  PSF,  and  bias  are  considered  maximum- 
likelihood  estimates.  Taking  the  conditional  expectation  of  Equation  (23)  results  in 

Q{pk,A,h,B)  = 

p°Jd  (to,  n;  Aold,p°kld,  hold )  In  [A  (to,  n)  pk  (to,  n)h(x-m,y  -  n)} 

k,x,y  ,m,n 

-  [A  (m,  n )  pk  (to,  n)h(x-m,y~  n)}  +  p°~d  ( x ,  y,  Bold )  In  [B  (x,  y)] 

- B(x,y )  (24) 

where 

pf d  (to,  n;  p? *,  Aold,  hold)  =  E 

and 

4d  (*>  26  B°ld)  =  E  [Qk  (*,  2/)  14  (*,  2/) ,  S°W]  •  (26) 

Equations  (25)  and  (26)  represent  the  expected  value  of  one  set  of  complete  data  given  the 
incomplete  data.  For  Poisson  random  variables,  these  expectations  turn  out  to  be  ratio  of  the 
data  times  one  expected  value  of  the  complete  data  divided  by  the  all  sets  of  expected  values  of 
the  complete  data  [11].  For  the  first  set  of  complete  data,  dk  (x.  y),  the  conditional  expectation 
is 


4  (x,  y\m,  n)  \dk  (x,y)  ,p°k 


old 


j^old  foold 


(25) 


4d  C m,n-,p°kld,A°ld,hold ) 


dk  (x,  y)  Aold  (to,  n)  pkd  (to,  n)  hold  (x  —  m,y  —  n) 
ikd  (x,y)  +  Bold  (x,y) 


(27) 


while  the  second  set  of  complete  data  concerning  the  pixel  bias  qk  (x.  y  ),  has  a  conditional 
expectation  equal  to 


\4d  yv\B°id) 


4  (x,y)  Bold  (x,y) 
ikd  (x,  y)  +  Bold  (x,  y ) 


(28) 


The  maximization  of  the  Q  function  for  all  parameter  unknowns  (target  amplitude,  target  pulse 
shape,  PSF,  and  pixel  bias)  is  theoretically  intractable  due  to  coupling.  It  is  required  to  break 
apart  the  Q  function  into  separate  components  such  that 


Q  —  Qp  +  Qh  +  Qa  +  Qb  (29) 

where  each  component  of  the  Q  function  can  be  maximized  independently.  Thus,  the  GEM 
algorithm  becomes 


QP(pTw\p°kd,Aold,hM) 

> 

( ^old\old  a  old  iold\ 

Qp  \Pk  \Pk  i  A  ,h  ) 

Qa  {Anew\p°kld,Aold:hold) 

> 

Qa  {Aold\p°kld,Aold:hold) 

Qh  (hnew\ptd1Aold1hold) 

> 

Qh  {hold\p°kld,Aold,hold) 

Qb  ( Bnew\Bold ) 

> 

Qb  ( Bold\Bold ) 

(30) 

which,  if  these  conditions  are  met,  ensures  that  the  likelihood  is  increased  with  each  iteration  [3] 
L  (pZw,Anew,  hnew,Bnew)  >  L  (plld,  Aold,  hold ,  Bold)  (31) 

resulting  in  a  GEM  sequence  converging  to  a  local  minimum. 
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Beginning  the  estimation  process  of  the  separate  <5  functions  starts  with  the  target  pulse 
shape,  Qp  which  is 


qp=  pT  (m> n;  Pfe 1 A°ldi  h°ld ) ln  [Pfc  (m> n)]  “  A  (m> n) 

kx,y,m,n 


K 


(m’n)  “  1 


,fc=i 


(32) 

where  a  pixel-dependent  Lagrange  multiplier,  A  (m,  n),  is  introduced  to  force  the  pulse  shape 
to  add  to  one  for  each  pixel.  This  constraint  is  necessary  to  decouple  the  pulse  shape  from  the 
target  amplitude  and  PSF.  Taking  the  derivative  of  Equation  (32)  with  respect  to  a  particular 
object  plane  point  and  range  sample,  setting  the  result  equal  to  zero,  dQp/dpka  ( m0 ,  n0)  =  0, 
and  solving  for  the  pulse  shape,  results  in 


PkT  (™0,  n0)  =  p°kl0d  ( m0 ,  n0)  — - ^ 

V  A  (m0,n0) 


(  Aold  (to0,  n0)\  {x,  y)  hold  ( x  -  m0,  y  -  n0) 

i°kld  {x,  V )  +  B°ld  (x,  y) 


(33) 


where 

A  ( rn0 ,  tl0 ) 

K 

=  Aold  (m0,n0)  ^2 

k=  1 

X  Y 


x—1 y— 1 


ifd  (x,  y )  +  Bold  (x,  y) 


and 


M  N 


i°kd  (x,y)  =  A°ld  (m’ n) pko  (m> n)  h°ld  ( x-m,y-n ) .  (35) 

m—  1  n—  1 

Equation  (33)  is  the  iterative  solution  for  the  pulse  shape  per  range  sample.  Next,  the  Q  function 
is  partitioned  into  its  target  amplitude  components 

M  N 

Qa=  ^2  {Pdd(m’n'’Pkd,Aold>hold)ln[A(m,n)]}-^'52A(rn,n)  (36) 

k,x,y,m,n  m—1  n= 1 


where 


J2J2h(x,y)  =  1 

x=l  y=  1 
K 

Y,pk  (m,n)  = 1 

fc= i 


(37) 

(38) 


have  been  utilized  to  decouple  the  pulse  shape  and  PSF  terms  from  the  target  amplitude.  Max¬ 
imizing  Equation  (36)  by  8Qa/9A  (toq,  n0)  =  0  and  solving  for  the  amplitude  term  results  in 
the  iterative  solution  for  the  target  amplitude  term 


A (m„, n„)  =  A-‘  (m„,  nc)±Pl“  (m„, »„)  £  £  dt  ~0, . 

(39) 

The  PSF  is  the  final  unknown  parameter  that  uses  the  first  set  of  complete  data,  dk  (x.  y).  The 
Q  function  for  the  PSF  is 


Qh  = 


E  Pf  (m’  A°ldi hM) ln  Mx-m,y-  n)] 

k,x,y,m,n 

-  [A  (to,  n)  pk  (to,  n)  h(x  —  m,y  —  n)\, 


(40) 
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which  still  has  the  target  amplitude  and  pulse  shape  terms.  Similar  to  [10],  a  change  of  variables 
is  required  to  remove  the  dependence  on  the  pulse  shape  and  to  allow  for  easier  differentiation. 
Utilizing  Yk=i  Pk  (m,  n)  =  1  and  setting  m!  =  x  —  m  and  n!  =  y  —  n,  Qh  then  becomes 


Qh 

V  {/if  [x-m',y-  n'-ptd ,  AM,  h°ld)  In  [h  (m',  n')]} 

k,x,y,m',n' 

—  ^  A(x  —  m' ,y  —  n')  h[rn! ,n').  (41) 

x,y,m',n' 

Setting  dQh/dh  ( m’0 ,  n^)  =  0  and  solving  for  the  PSF  produces  the  iterative  solution 

hnew  ,,  /  X  _  h  old  ,,  /  \  V-  4  (a;,  y)A{x-m'0,y-  n'0)  p°kld  (x-m'0,y-  n'0) 

ii  yiii0)  fi0)  ii  yiri0)  n0j  y  ^  ^ 

k,x,y  (ioM  (a.)  +  Bold  (X,  y))  Y  Y  A{X~  TOg,  y  -  <) 

£C=1  y— 1 

(42) 

Usually,  the  target  amplitude  term  in  the  denominator  would  be  an  issue  because  it  is  consid¬ 
ered  the  new  estimate.  However,  Equation  (39)  is  the  new  estimate  and  can  replace  the  target 
amplitude  in  the  denominator  resulting  in  a  consistent  solution  for  the  PSF.  Finally,  the  pixel 
bias  must  be  estimated.  In  order  to  estimate  the  pixel  bias,  the  second  set  of  complete  data, 
qk  (x,  y),  is  utilized.  The  Q  function  for  the  pixel  bias  is 

=  Ilf;  -  '  ■"  (B  (*.  y))-B  (X, ,).  (43) 

,  il  y)  +  B  y ) 

k=l  x=l  y—1  K  v  '  ' 

Setting  dQs/dB  ( x0 ,  y0)  =  0  and  solving  for  the  pixel  bias  results  in  an  iterative  solution 

Bnew  (Xo,  Vo)  =  B°ld  (Xo,  Vo)  J2  ,.old(  - yT  (44) 

£ri  [}k  ( xo ,  Vo)  +  Bold  (x0,  y0)) 


corresponding  to  the  pixel  bias  solution. 

After  a  pre-determined  number  of  iterations  on  Equations  (33),  (39),  (42),  and  (44),  range 
estimate  updates  for  each  pixel  are  generated  by  using  a  normalized  correlation  method  between 
a  reference  waveform  at  sub-sample  ranges  and  the  the  GEM  estimate  for  the  pulse  shape,  p%ew . 
The  range-dependent  reference  waveform  that  results  in  the  highest  correlation  is  chosen  and 
the  corresponding  range  is  the  new  range  estimate  for  that  pixel.  The  new  range  estimate  is  fed 
back  into  the  pulse-shape  to  generate  a  new  p";!d  followed  by  another  set  of  GEM  iterations.  The 
process  (GEM  iterations  followed  by  range  updates)  repeats  with  the  new  range  estimates  used 
in  calculating  p^ld  using  Equation  (33)  and  ceases  when  the  mean  square  error  (MSE)  between 
the  data  and  non-noisy  range  images  reaches  the  stopping  criteria.  All  previous  amplitude, 
PSF,  and  pixel  bias  estimates  carry  over  from  one  range  update  to  the  next.  More  specifically, 
iterations  cease  when  the  MSE  is  lower  than  the  average  data  variance  or 

K  X  Y  K  X  Y 

EEE (dik (x’y)~- Tkst (x,y)~ Bnew Ou y)) 2  <  E]  Vk (x’ y)  (45) 

k—1  x=l  y—1  k=l  x=l  y—1 

with 

M  N 

Iekst  0 * ,  !/)=EE  Anew  (m’ n)  PnkeW  (m’ n)  hnew  ( x-m,y-n )  (46) 

m—  1  n=  1 
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and 


/ 


2 


Vk  {x,y)  =  Y 
j=  i 


djk  (x,  y)  - 


E  dh  k{x,y) 

32  =  1 _ 

J 


\ 


(47) 


where  the  extents  of  j2  are  the  same  as  j  and  14  is  an  average  variance  for  the  volume  elements 
(voxels)  across  many  data  cubes.  In  the  experimental  data,  a  specific  distribution  for  the  variance 
is  not  chosen  in  order  to  account  for  all  noise  sources.  For  the  simulation  data,  the  data  variance 
is  replaced  by  the  average  variance  of  the  Poisson  distribution  due  to  it  being  a  known  and 
controlled  noise  source. 

Once  the  GEM  estimates  the  pulse-shape,  the  range  for  each  pixel  must  be  extracted.  If 
simple  peak  detection  is  used,  targets  between  the  samples  would  have  the  wrong  range  es¬ 
timate.  In  order  to  mitigate  inter-sample  targets,  scaling,  and  waveform  truncation  issues, 
sub-sample  ranging  is  performed  on  the  pulse-shape  by  using  a  normalized  correlation  method 
based  on  the  Pearson  product-moment  correlation  coefficient.  Using  this  coefficient  forces  each 
pixel’s  waveform  to  be  zero  mean  and  unit  standard  deviation.  The  normalized  cross-correlation 
method  is  constructed  as  follows:  The  range  vector  of  samples  within  a  cube  R  (t)  is  represented 
by 

K- 1 

K(l)  =  £  (•^minH^inc  (f))  (48) 

t= 0 


where  K  is  total  number  of  samples,  zm[n  is  the  range  of  the  first  sample,  and  z-mc  is  the  range 
increment  per  sample.  Another  range  vector,  Kr  ( t ')  is  constructed  with  the  same  maximum 
and  minimum  extents  as  R  (t),  but  with  a  smaller  range  increment  per  sample  or 


K'-l 

Kr  (t)  =  Y  (^min+Zf  (t)) 

t—0 


(49) 


where  K'  is  the  number  of  samples  in  Kr  and  Zf  is  the  range  increment.  Since  the  extents  of 
Kr  match  R(p),  K'  >  K  and  Zf  <  z\nc.  A  2D  reference  Gaussian  waveform  matrix  is  used 
with  the  Kr  vector  as  the  reference  target  location  or 


r k  (p)  =  exp 


-(4 


2 Kr  ( p )  /c)2 

2<4 


(50) 


where  tk  is  the  time  vector  and  k  £  [1,  K] 

S2  ( P ) 

where  of  and  fk  are  the  variance  and  average  of  rk  in  the  time  dimension.  Considering  the 
range  estimate  for  the  (to,  n)th  pixel,  the  zero  mean  and  unit  variance  version  of  the  estimated 
pulse-shape,  pk  (to,  n )  is 


The  zero  mean  and  unit  variance  version  of  rk  is 

rk  (p)  ~  fk  (p) 

V2r(p) 


Si  (to,  n)  = 


pk  (to,  n)  - 


K 

Y  Pk(rn,n) 
k= 1 

K' 


(m,  n) 


(52) 


where  <t?“  is  the  variance  of  pk  (to,  n)  in  the  time  dimension.  With  Si  and  S2  determined,  the 
normalized  cross  correlation  denoted  by  *  is  performed  by 

S2*Si 


CKr  ( V )  = 


K' 


(53) 
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and  finding  the  range  estimate  is  accomplished  by  peak  detection  on  Ckt  or 

R(m,n)  =  arg  max  Cnr  ( p )  (54) 

2f(p) 

In  summary,  the  pulse-shape  estimation  algorithm  steps  are: 

1 .  Initialize  PSF,  amplitude,  and  pixel  bias 

2.  Determine  initial  ranges  and  input  into  pulse-shape 

3.  Perform  GEM  iterations  using  Equations  (33),  (39),  (42),  and  (44) 

4.  Use  normalized  cross-correlation  to  find  new  range  estimates  with  Equation  (54) 

5.  Generate  new  pulse-shapes  based  on  new  ranges 

6.  Determine  MSE  and  compare  to  stopping  criteria 

7.  Repeat  Steps  3  through  6  until  stopping  criteria  violated 

8.  Range  estimates  taken  from  last  execution  of  Step  4 


2.4  Range  estimation  using  object  recovery  via  the  GEM  algorithm 

When  multiple  cubes  are  available  and  properly  registered  spatially  and  temporally,  another 
method  to  perform  range  estimation  is  to  relax  the  constraint  on  the  pulse-shape  and  assume 
just  an  object  in  the  data  model.  This  change  mitigates  the  issue  in  the  hardware  data  where  the 
pulse-shape  is  vaguely  known.  Therefore,  estimation  is  performed  on  O/-  rather  than  on  pk  from 
Equation  (8).  The  problem  setup  is  similar  to  the  pulse-shape  estimation  (now  with  more  than 
one  cube)  by  calling  the  original  data,  djk{x ,  y),  the  incomplete  data  and  specifying 

M  N 

djk  {x,  y)  =  EE  djk  (x,  t/|rn.,  fi)  T  Qjk  ( x ,  y)  (55) 

m=  1  n=  1 


where  two  new  variables,  djk(x,  y\m ,  n )  and  qjk  ( x ,  y),  are  defined  and  called  complete  data. 
This  formulation  provides  two  sets  of  complete  data  that  account  for  the  image  formation  and 
pixel  bias  respectively.  The  same  PSF  can  be  assumed  for  adjacent  collections  due  to  a  typical 
data  collection  scenario  where  environments  shouldn’t  be  changing  rapidly  (ignore  j).  Thus, 
the  expected  values  of  the  complete  data  sets  are  given  by 


E 


djk  ( x,y\m,n ) 


Ok  ( m ,  n)  h  [x  —  m,y  —  n) 


(56) 


and 

E  [qjk  (. x ,  y)]  =  B  ( x ,  y)  (57) 

where  B  (x,  y)  is  the  constant  pixel  bias.  The  expected  value  of  the  incomplete  data  is  thus 

E  [d]k  (x,  y)]  =  ik  (x,  y)  +  B  (x,  y) .  (58) 

The  PMF  for  the  photon  noise  is 


P 


m.n)\  = 


(dk  {x,y  | 

[ok  (■ to ,  n)  h(x  —  m,y  —  n)]^fc(x,y lm’")  exp  {— Ok  (m,  n)  h  (x  —  m,y  —  n)} 

djk  (x,  y\m,  n)\ 


(59) 
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while  the  pixel  bias  PMF  is 


P(Qjk  (x,y)) 


B(x,yyjk(x’v)e~B^ 

Qjk  (x,y)\ 


(60) 


Assuming  statistical  independence  between  all  the  pixels  and  between  the  photon  noise  and 
pixel  bias  noise,  the  complete  data  log-likelihood  is  then 


or 


Lcd  (ofc,  h,  B)  =  In 


n 

j,k,x,y,m,n 


(djk  (x,  y\m ,  n)j  P  (qjk  (x,  y)) 


(61) 


LCd  (ok,  h,  B)  =  ^2  djk(x,y\m,n)ln[ok{m,n)h(x-m,y-n)] 

j,k,x,y,m,n 

-  [ok  (to,  n)h(x  -  m,y  -  n)]  +  qjk  {x,  y)  In  [B  (x,  y)]  -  B  (x,  y) .  (62) 


Referring  to  [3],  the  Q  function  becomes  the  expected  value  of  the  complete  data  log-likelihood 
function  with  respect  to  the  incomplete  data  and  old  parameter  estimates 

Q  0 ok ,  h,B)  =  E  [LCd  (ok,  h,  B)  \ djk  {x,  y) ,  of,  hold ,  B°ld]  .  (63) 

Taking  the  conditional  expectation  from  Equation  (63)  results  in 

Q  (ofc,  h,  B)  =  ^2  (TOi  n\of,hold)  In  [ok(m,n)h  (x  -m,y-n)] 

j,k,x,y  ,m,n 

-  [ok  (m,  n)h(x-m,y-  n)}  +  y°~d  (at,  y\ Bold )  In  [B  (x,  y)]  -  B  (x,  y )  (64) 


where 


y°~d  (m,  n;  okld,  hold )  =  E  djk  (x,  y\m,  n)  \djk  ( x ,  y) ,  of,  hold 

djk  (x,  y)  of  (to,  n)  hold  (x  —  m,y  —  n) 
if  (x,  y)  +  Bold  (x,  y) 


and 


(65) 


df  {x,y\Bold)  = 


E  [qjk  (x,  y)  |  djk  {x,  y) ,  Bold\ 
djk  (x,y)  Bold  (x,y) 
if  (x,  y)  +  Bold  (x,  y) 


(66) 


Similar  to  the  pulse-shape  estimation,  the  maximization  of  the  Q  function  for  all  parameter 
unknowns  (object,  PSF,  and  pixel  bias)  is  theoretically  intractable  due  to  coupling.  It  is  required 
to  break  apart  the  Q  function  into  separate  components  such  that 


Q  —  Qo  +  Qh  +  Qb 


(67) 


where  each  component  of  the  Q  function  can  be  maximized  independently.  Thus,  the  GEM 
algorithm  becomes 


Qo  (oTw\of,h°ld) 
Qh  (hnew\of,hold) 
Qb  (Bnew\Bold) 


>  Q0(of\of,h°ld) 

>  Qh(hold\af,hold ) 

>  Qb  (Bold\Bold) 


(68) 
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(69) 


ensuring  that  the  likelihood  is  increased  with  each  iteration  [3] 

L  {olew,hnew,  Bnew)  >  L  ( o°kld7  hold:  Bold ) 

resulting  in  a  GEM  sequence  converging  to  a  local  minimum.  The  procedure  to  find  the  maxima 
of  the  Q  functions  is  the  same  as  in  pulse-shape  estimation.  First,  the  object  solution  is  found 
by  specifying 

Qo  =  5Z  H°Jd  (m,  n]0%ld,hold)  In  \ok  (m,n)\  -  ok  (m,n)  h  (x  -  m,  y  -  n) .  (70) 

j,k,x,y,m,n 

In  order  to  maximize  Qa,  the  derivative  of  Equation  (70)  with  respect  to  a  particular  object 
plane  point  and  range  sample  is  set  equal  to  zero,  dQp/doka  ( m0,n0 )  =  0.  Solving  for  the 
object  results  in 


„new  \  _  °°kt  (m°’  n° )  V'  djko  {%,  y)  hold  ( X  —  Too,  y  —  Tl0) 

°i.  ( m0 ,  n0)  —  T  /  y  y  y ' 


J 

j=i  *=i  y=t 

with  J  as  the  number  of  cubes  and  utilizing 


i^d  (x,  y)  +  Bold  (x,  y) 


(71) 


Y 

(x--  y^> =  1 

x—l  y—1 


(72) 


and  where 

M  N 

i°kd  {x,  v)  =  '52'52  °kld  (m> n)  h°ld  (x~m,y-  n )•  (73) 

m= 1 n— 1 

Equation  (71)  is  the  iterative  solution  for  the  pulse  shape  per  range  sample.  The  PSF  is  the  other 
unknown  parameter  that  uses  the  first  set  of  complete  data,  d:jk  (x.  y).  The  Q  function  for  the 
PSF  is 


Qh  =  y°^ld  (m,n-olld,  hold)  In  [h(x-m,y -n)}- ok  (m,n)h(x-m1y-n). 

j,k,x,y,m,n 

(74) 

Similar  to  [10],  a  change  of  variables  is  required  to  remove  the  dependence  on  the  pulse  shape 
and  to  allow  for  easier  differentiation.  Setting  m!  =  x  —  m  and  n1  =  y  —  n,  Qh  then  becomes 


Qh  =  5Z  d  ( x  ~  m' j  V  ~  n'i  °kdi  hold)  In  [h  (■ m n')]—Ok  {x  —  m' ,  y  —  n')h  (to7,  n') 

j  ,k  ,x  ,y  ,m'  ,n' 

(75) 

Setting  dQh/dh  (m'0,  n'Q)  =  0  and  solving  for  the  PSF  produces  the  solution 


j-new 


{m'0,n'0) 


h°ld  (m'0,  n'0) 


J 


E  Ofc  (x  —  m'0,  y  —  n'0) 

k,x,y 


E 


djk  (x,  y)  o°kld  (x  -m'0,y-  n'Q) 
ifd  (x,  y)  +  Bold  (x,  y) 


(76) 

The  object  term  in  the  denominator  is  the  new  estimate  from  Equation  (71).  Since,  there  are 
phase  abberations  across  the  aperture  and  the  PSF  needs  to  be  constrained,  phase  retrieval  is 
performed  on  Equation  (76)  by  the  Gerchberg-Saxton  algorithm  [12].  In  the  pulse-shape  es¬ 
timation,  it  was  the  object  (i.e.  pulse-shape)  that  was  constrained  making  the  phase  retrieval 
unnecessary.  Finally,  the  pixel  bias  must  be  estimated.  In  order  to  estimate  the  pixel  bias,  the 
second  set  of  complete  data,  qk  (x.  y),  is  utilized.  The  Q  function  for  the  pixel  bias  is 


Qb 


K  X  Y 


EEE 


djk  (x1y)Bold  (x,y) 
i%ld  (x,  y)  +  Bold  (x,  y) 


ln(S  (x,y)) 


B  (x,  y). 


(77) 
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Table  1.  3D  FLASH  LADAR  parameters 


Parameter 

Value 

Detector  Array 

128  x  128 

Aperture  Diameter  ( D ) 

2  mm 

Mean  Wavelength 

1.55  pm 

Focal  Length 

0.30  m 

Target  Range 

5.21  m 

Transmit  Energy 

10  mJ 

Pulse  Standard  Deviation  ( aw ) 

3  ns 

Beam  Divergence 

0.009  radians 

Detector  Size 

10  pm 

Detector  Spacing 

100  pm 

Detector  Array  Fill  Factor 

100% 

Detector  Bandwidth 

0.5  pm 

Target  Reflectivity 

10% 

Detector  Array  Size 

30  x  30  pixels 

Solar  Irradiance 

10  Watts/m2/ pm 

D/r0  Seeing  Condition 

1.43 

Frame  Rate 

30  Hz 

Time  Samples  (per  collect) 

20 

Sample  Period  (within  a  collect) 

1.876  ns 

Setting  dQs/dB  (x0,  y0 )  =  0  and  solving  for  the  pixel  bias  results  in 


■Qnew 


(x0,y0) 


Bold  (x0,y0) 
JI< 


djk  {x0,  yo ) 


{x0,  y0)  +  Bold  (x0,  ya) 


(78) 


GEM  iterations  continue  and  cease  when  the  mean-square  error  (MSE)  violates  the  stopping 
criteria.  Once  the  stopping  criteria  is  reached,  range  estimates  are  determined  by  using  the 
normalized  cross-correlation  method  on  the  object  estimate  as  described  in  the  previous  section. 
In  summary,  the  object  estimation  steps  are: 

1 .  Initialize  object,  PSF,  and  pixel  bias 

2.  Perform  one  GEM  iteration  using  Equations  (71),  (76),  and  (78) 

3.  Determine  MSE  and  compare  to  stopping  criteria 

4.  Repeat  Steps  2  and  3  until  stopping  criteria  reached 

5.  Use  normalized  cross-correlation  to  find  new  range  estimates  with  Equation  (54) 


3  SIMULATION 

In  order  to  verify  the  theory,  a  simulation  scenario  was  developed  whereby  targets  are  interro¬ 
gated  by  a  3D  FLASH  LADAR  defined  by  the  parameters  from  Table  1 .  The  goal  is  to  perform 
range  estimation  given  the  noisy,  blurry  data  observations.  Results  show  range  estimation  im¬ 
provement  by  performing  object  recovery  either  via  a  Wiener  filter  method  or  GEM  algorithms 
as  outlined  in  Sections  2.3  and  2.4.  Our  previous  research  has  taken  the  approach  to  use  a 
Wiener  filter  on  each  individual  range  slice  and  then  use  a  pixel-based  ranging  method  on  the 
resulting  “deblurred”  data  cube  [5].  Performance  will  illustrate  that  the  GEM  algorithms  pro¬ 
vide  increased  error  performance  over  the  Wiener  filter  while,  at  the  same  time,  being  more 
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Fig.  4.  True  ranges  for  simulation  targets:  (a)  three  bars,  (b)  Many  bars,  (c)  Various  blocks,  (d)  Cylinder, 
(e)  Slanted  boards,  and  (f)  Connected  blocks.  The  target  names  in  this  caption  correspond  to  the  targets  in 
Table  2.  The  three  bar  target  is  also  the  experimental  data  target.  Other  targets  illustrate  the  robustness  of 
the  estimation  algorithms. 
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Table  2.  Range  estimation  results  for  simulation  data 


Data  set 

Metric 

OD 

WF 

GEMp 

GEM0 

Three  bars 

RMSE  (m) 
Corr 

0.402 

0.767 

0.346 

0.830 

0.163 

0.963 

0.100 

0.984 

Many  bars 

RMSE  (m) 
Corr 

0.596 

0.687 

0.561 

0.664 

0.346 

0.786 

0.365 

0.794 

Slanted  boards 

RMSE  (m) 
Corr 

0.225 

0.945 

0.171 

0.971 

0.161 

0.967 

0.131 

0.983 

Cylinder 

RMSE  (m) 
Corr 

0.184 

0.877 

0.153 

0.925 

0.160 

0.945 

0.153 

0.962 

Various  blocks 

RMSE  (m) 
Corr 

0.473 

0.595 

0.209 

0.931 

0.344 

0.725 

0.175 

0.955 

Connected  blocks 

RMSE  (m) 
Corr 

0.208 

0.853 

0.133 

0.955 

0.158 

0.918 

0.112 

0.970 

robust.  Again,  the  GEM  algorithms  are  more  robust  in  that  they  do  not  need  to  know  the  point 
spread  function,  unlike  the  Wiener  filter  technique. 

Using  a  Gaussian  transmitted  pulse,  a  3D  FLASH  LADAR  imaging  scenario  is  developed 
in  simulation  using  various  geometrical  shapes  as  targets  shown  in  Fig.  4(a)-(f).  One  important 
clarification  on  the  receiver  optics  is  that  the  detector  array  has  an  effective  fill  factor  of  100% 
by  placing  a  micro-lens  array  in  front  of  the  pixels  to  focus  the  light  onto  the  pixel.  Also,  the 
data  includes  effects  from  an  average  atmospheric  turbulence  to  enable  blind  deconvolution. 
Range  estimates  are  also  determined  without  processing  to  enable  further  comparison  between 
no  processing  and  object  recovery  attempts.  Results  for  all  the  targets  and  methods  with  error 
metrics  are  summarized  in  Table  2. 

Table  2  clarifications:  “RMSE”  is  root  mean  square  error  (RMSE)  in  meters  between  the 
true  ranges  and  estimated  ranges,  “Corr”  is  an  image  quality  metric  referring  to  the  correlation 
between  the  true  range  image  and  estimated  range  image  signifying  linear  relationship  strength 
(not  to  be  confused  with  the  normalized  correlation  ranging  method),  “OD”  refers  to  the  orig¬ 
inal  data  (OD)  with  no  deblurring  and  ranges  estimated  by  the  normalized  cross-correlation 
method,  “WF”  relates  to  range  estimation  using  a  Wiener  Filter  technique  with  normalized 
cross-correlation  [5],  “GEMp”  is  the  pulse-shape  estimation  GEM  algorithm,  and  “GEM0”  is 
the  object  estimation  GEM  algorithm. 

The  targets  of  primary  interest  are  the  three  bar  target  and  the  multiple  bar  target  because 
the  three  bar  target  is  also  the  experimental  target  and  the  multiple  bar  target  is  most  sensitive 
to  spatial  blurring  of  all  the  targets.  The  bar  targets  are  constructed  in  simulation  consisting  of 
two  flat,  perpendicular  optically  rough  surfaces  at  different  ranges.  Referring  to  Figs.  4(a)  and 
(b),  the  first  surface  in  range  has  rectangular  cut-out  shapes  while  the  second  surface  contains 
no  cutouts.  This  type  of  target  was  chosen  to  highlight  not  only  the  coupling/blurring  effects 
of  the  pixels  along  the  edges  of  the  rectangles,  but  also  the  decoupling  and  ranging  capability 
of  the  GEM  algorithm.  The  other  targets  are  built  in  similar  manner.  Bar  target  shapes  were 
used  because  the  distances  and  shape  dimensions  can  be  physically  measured  in  a  laboratory 
environment  to  show  range  estimation  improvement. 

Table  2  and  the  range  images  from  Figs.  4  and  5  show  the  negative  effects  of  the  blurring 
on  range  estimation  juxtaposed  with  the  positive  effects  from  attempting  recover  the  original 
object  through  Wiener  filtering  or  the  GEM  algorithms.  Implicit  in  the  results  is  the  ability 
to  accurately  estimate  the  pixel  bias.  Without  it,  the  object  model  falls  apart  and  range  error 
becomes  extremely  large.  Through  simulation,  the  model  and  object  recovery  attempts  have 
been  verified.  The  final  step  is  to  use  experimental  data  to  validate  simulation  results. 
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Fig.  5.  Estimated  ranges  for  simulation  targets:  (a)  No  processing  -  three  bars,  (b)  GEM0  processing  - 
three  bars,  (c)  No  processing  -  Many  bars,  and  (d)  GEM0  processing  -  Many  bars.  Utilizing  the  GEM0 
algorithm,  simulation  results  show  the  image  quality  improvement  and  improved  range  estimation  ( RMSE 
improves  75%  for  3  bar  target). 
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Fig.  6.  (a)  Using  the  data  from  Fig.  5(c)-(d),  investigating  pixel  (8,23)  shows  the  estimated  waveform 
(object  plus  pixel  bias)  closely  matching  the  true  waveform  while  the  detected  waveform  does  not.  The 
estimated  range  is  6.73  m  while  the  true  range  is  6.7  m.  The  algorithm  also  implicitly  estimates  the  pixel 
bias  term  accurately,  (b)  Again,  using  the  data  from  Fig.  5(c)-(d),  investigating  pixel  (16,17)  shows  the 
estimated  waveform  improving  upon  the  detected  waveform,  but  not  able  to  match  the  true  waveform  as 
well  as  the  previous  pixel.  The  estimated  range  is  6. 1 1  m  while  the  true  range  is  6.7  m.  Incorrect  range 
estimation  after  the  GEM0  algorithm  relates  to  blurring  severity  (edges  of  cut-outs  in  first  surface)  and/or 
a  particularly  noisy  realization  from  the  Poisson  distribution. 
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4  EXPERIMENTAL 


With  the  simulation  giving  credence  to  the  theory,  experimental  results  are  presented  to  fully 
validate  the  data  model  and  the  claim  of  improved  range  estimation  in  3D  FLASH  LADAR 
via  object  recovery.  Using  the  3  bar  target  template,  a  laboratory  experiment  was  conducted 
using  3D  FLASH  LADAR  hardware  from  AFRL  consistent  with  parameters  in  Table  1 .  As  in 
simulation,  experimental  results  show  range  estimation  improvement  after  applying  the  object 
recovery  techniques.  However,  several  modifications  to  the  camera  and  raw  data  were  necessary 
to  enable  a  proper  experiment. 


4.1  Spatial  aliasing 

Due  to  limits  in  current  detector  technology  requiring  a  large  footprint  for  the  electronics  behind 
each  pixel,  the  receiver  optics  are  spatially  under-sampled  which  needs  to  be  mitigated  in  order 
for  the  received  data  to  be  unaliased.  This  determination  comes  from  Nyquist  sampling  theory 
in  which  the  sampling  rate  must  be  at  least  twice  the  highest  frequency  content  in  the  signal. 
The  optics  are  a  natural  low-pass  filter  with  the  highest  frequency  called  the  cut-off  frequency. 
For  incoherent  imaging,  the  cut-off  frequency  is  [6] 


fo 


D 

Xzi 


(79) 


where  D  is  the  aperture  (exit  pupil)  diameter,  A  is  the  light  wavelength,  and  z,  is  the  image 
distance.  Therefore,  the  focal  plane  must  sample  at  twice  this  spatial  frequency  or  2D/Xzi. 
The  typical  apertures  for  this  camera  are  in  the  centimeters.  For  example,  an  aperture  of  10  cm 
would  equate  to  a  spatial  frequency  sampling  requirement  at  4.3xl05  cycles  per  meter.  At  100 
/im  spacing,  the  detector  array  does  not  meet  this  requirement.  If  the  aperture  is  reduced  to  2 
mm,  then  the  spatial  frequency  sampling  requirement  is  now  at  8.6xl03  cycles  per  meter  which 
the  detector  array  can  meet.  However,  the  aperture  reduction  comes  at  the  expense  of  reduced 
light  gathering  and  shortened  range  in  which  the  LADAR  can  be  operated.  Thus,  the  target 
range  is  placed  at  5.21  meters  (near  the  minimum  ranging  distance  of  the  sensor)  to  obtain  high 
enough  signal  to  noise  ratio  (SNR)  in  the  collected  data. 


4.2  Data  pre-processing 

The  data  observations  from  the  3D  FLASH  LADAR  hardware  need  pre-processing  steps  to  be 
suitable  for  insertion  into  the  Wiener  filter  and  GEM  algorithms.  In  simulation,  the  noisy  and 
blurry  data  is  well-controlled  and  therefore,  well-behaved.  While  the  experimental  3D  FLASH 
LADAR  data  exhibits  expected  pixel  waveform  shapes  (i.e.  Gaussian-like)  and  spatial  blur,  the 
data  is  ill-behaved  to  a  degree  due  to  inherent  features  of  the  hardware  performance. 

Referring  to  [13]  and  [14],  the  experimental  hardware  experiences  a  gain  phenomenon 
whereby  a  pixel’s  gain  drops  when  laser  energy  is  incident  upon  a  large  area  of  another  part 
of  the  detector  array.  With  the  3  bar  target,  the  laser  energy  is  incident  front  surface  first  which 
causes  second  surface  pixels  to  experience  a  gain  drop.  Figure  7(b)  shows  the  gain  drop  for  a 
second  surface  pixel.  The  method  for  correcting  the  gain  is  to  calculate  an  average  gain  profile 
by  looking  at  background  pixels  (i.e.  returned  laser  energy  not  incident  on  these  pixels). 

Assuming  the  system  noise  follows  the  Poisson  distribution  and  the  gain  is  constant  between 
pixels,  the  data  model  for  an  arbitrary  pixel  is 

d(t)  =  G  (t)  [Js  (f)  +  IB  (t)]  (80) 

where  G  ( t )  is  the  unitless,  time-varying  gain.  Is  (t)  is  the  laser  signal  in  units  of  photons,  and 
Ib  ( t )  is  the  background  signal.  A  new  variable  d(t)  is  determined  by 

d{t)  =  $^-  (81) 

(i) 
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(b) 


Fig.  7.  (a)  Gain  profile  correction  resulting  from  executing  Equation  (84).  By  looking  at  background 
pixels,  the  hardware  gain  dip  is  clearly  evident  at  the  first  surface  (near  range  sample  five)  and  the  second 
surface  (near  range  sample  nine).  The  first  surface  gain  drop  is  larger  than  the  second  surface  gain  drop 
due  to  the  larger  number  of  pixels  illuminated  (i.e.  larger  surface  area).  Amount  of  gain  drop  is  propor¬ 
tional  to  received  intensity  level  and  quantity  of  pixels  illuminated,  (b)  Investigating  Pixel(  19,32)  from 
experimental  three  bar  target,  the  pixel  waveform  benefits  from  the  gain  variation  correction  by  removing 
the  gain  drop  near  range  sample  four.  After  correction,  the  pixel  waveform  looks  more  like  the  intended 
pulse  model,  but  with  unwanted  noise  artifacts. 
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where  is  if)  =  G  (t)  I b  (t)  and  is  a  known  average  background  signal.  It  is  separately  calcu¬ 
lated  in  the  laboratory  by  averaging  the  detected  background  signal  for  selected  voxels  across 
many  data  cubes.  Looking  at  the  background  pixels  only,  d  (t)  is 


d  (t)  = 


G  (t)  IB  (t)  G  (t)  IB  (t)  I  sit) 


%b  it)  G  (t)  Is  it)  Ib  it) 
Taking  the  statistical  variance  results  in 

Ib  it)  Ib  (t)s 


(82) 
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Applying  this  result  and  using  a  sample  variance  of  d  in  place  of  the  statistical  variance  (s2  — > 
var  (d  (f  ),  the  gain  is  determined  by 

Git)  =  T~m  =  iB  (t)s2-  (84) 

Ib  it) 

and  can  be  seen  in  Fig.  7(a).  This  gain  profile  is  used  on  each  of  the  pixels  waveforms  to  correct 
for  the  hardware  deficiencies  and  to  more  closely  match  the  model.  For  example.  Fig.  7(b) 
shows  the  benefits  of  the  gain  correction  for  one  second  surface  pixel.  Also  observed  in  the 
previous  work,  a  side  benefit  of  gain  correction  in  both  first  and  surface  pixels  is  the  waveform 
becomes  more  symmetrical.  The  emitted  laser  pulse  shape  is  a  hybrid  of  a  Gaussian  or  negative 
parabolic  shape  with  some  asymmetry.  Gain  correction  takes  out  some  of  the  asymmetry. 

The  3D  FLASH  LADAR  is  also  not  a  photon-counting  device  where  one  digital  count  equals 
one  photon.  The  receiver  optics  use  Avalanche  Photo  Diodes  (APD)  where  one  photon  equals 
many  detected  counts.  Consequently,  intensity  scaling  must  be  performed  to  condition  the  data 
to  be  consistent  with  the  Poisson  distribution.  The  conditioning  is  performed  by  using  the 
statistics  of  the  light  and  the  detected  mean  and  variance  of  the  data.  The  detected  mean  of  the 
data  is  qK  where  q  is  a  scaling  factor  with  units  of  photons  per  detected  counts  and  K  is  the  true 
mean  in  units  of  photons.  Since  incoherent  imaging  is  assumed,  the  detected  variance  becomes 

gV  =  q2  K  (85) 


noting  that  the  mean  and  variance  of  the  Poisson  distribution  are  the  same.  The  data  is  scaled 
by  solving  for  q  and  then  converting  the  detected  counts  to  photons  by 


(86) 


where  dph  is  the  data  in  units  of  photons  and  d.,ic  is  the  data  in  units  of  detected  counts. 


4.3  Experimental  PSF 

The  Wiener  filter  is  used  to  provide  a  comparison  to  the  GEM  algorithm  [5].  In  order  to  imple¬ 
ment  the  Wiener  filter,  the  PSF  must  be  known.  Since  the  derivative  of  a  system  step  response  is 
the  system  impulse  response,  the  PSF  is  determined  by  taking  the  derivative  of  a  experimental 
step  target.  Figure  8(a)  shows  a  range  image  of  the  step  target  collected  with  the  same  hardware 
as  the  bar  target  data.  Although,  the  entire  range  image  does  not  meet  the  requirements  of  being 
a  step  target  due  to  the  non-uniform  intensity  on  the  left-hand-side  (LHS).  Therefore,  a  sym¬ 
metric  impulse  response  was  assumed  and  the  right-hand-side  (RHS)  of  the  impulse  response 
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Fig.  8.  (a)  One  range  image  of  the  step  target  data  cube.  Although  the  board  edge  is  clearly  visible,  the 
variable  intensity  across  it  causes  an  issue  with  the  impulse  response  calculation.  The  step  response  defi¬ 
nition  requires  a  constant  amplitude  at  all  spatial  positions.  The  target  board  portion  of  the  step  response 
does  not  meet  this  requirement,  but  the  non-target  area  (right-hand-side)  does  exhibit  a  constant  ampli¬ 
tude.  The  portion  of  the  step  response  function  where  it  turns  off  is  this  non-target  area.  Performing  the 
step  response  derivative  only  on  this  non-target  area  solves  the  problem  of  variable  target  board  amplitude, 
(b)  ID  cut-out  of  the  resulting  PSF.  Assuming  circular  symmetry,  an  outer  product  operation  is  used  to 
find  the  corresponding  2D  PSF  function,  (c)  Optical  transfer  function  (OTF).  The  OTF  is  found  by  taking 
Fourier  Transform  of  the  experimental  PSF  [6],  (d)  ID  cut-out  (zero  spatial  frequency)  of  the  OTF.  The 
profile  shows  nearly  diffraction-limited  optics  with  a  cut-off  frequency  at  4050  cycles  per  meter. 


was  copied  and  flipped  over  to  use  as  the  LHS.  Figure  8(b)  exhibits  the  resulting  profile  with  an 
outer  product  operation  producing  the  two-dimensional  PSF.  Phase  retrieval  is  then  performed 
via  the  Gerchberg-Saxton  algorithm  to  arrive  at  the  PSF  used  by  the  Wiener  filter  [12].  This 
requirement  to  know  the  PSF  is  a  shortcoming  of  the  Wiener  filter  algorithm.  Figures  8(c)-(d) 
show  the  optical  transfer  function  (OTF)  where  the  optics  exhibit  a  nearly  diffraction-limited 
performance. 

4.4  Results 

Table  3  and  Fig.  9  illustrate  the  range  estimation  benefits  of  object  retrieval.  The  pulse-shape 
and  object  estimation  give  an  RMSE  improvement  of  25%  and  34%  respectively  over  the  orig¬ 
inal  data.  Additionally,  the  pulse-shape  and  object  estimation  give  an  RMSE  improvement  of 
7%  and  18%  respectively  over  the  Wiener  filter  algorithm.  Figure  9(c)  shows  the  image  quality 
improvement  over  the  original  data  range  image  in  Fig.  9(b).  Pixel  waveforms  provide  addi¬ 
tional  information  on  the  object  recovery  abilities.  Figure  9(d)  demonstrates  this  ability  on  a 
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Table  3.  Range  estimation  results  for  experimental  data 


Data  set 

Metric 

OD 

WF 

GEMp 

GEM0 

3  bars 

RMSE 

Corr 

0.301 

0.818 

0.243 

0.883 

0.226 

0.900 

0.198 

0.924 

(b) 


pixel 

(c) 


Fig.  9.  Experimental  target:  (a)  True  ranges  with  first  surface  at  5.21  m  and  second  surface  at  6.43  m  with 
1 .22  m  of  separation  in  between  surfaces,  (b)  Ranges  using  normalized  cross-correlation  without  using 
estimation,  (c)  Estimated  ranges  using  GEM0  algorithm  followed  by  normalized  cross-correlation,  (d) 
Considering  pixel  (32.18),  its  estimated  waveform  (object  plus  pixel  bias)  shows  similar  results  from  the 
simulated  data.  The  estimated  waveform  more  closely  resembles  the  true  waveform  with  the  range  close 
to  range  sample  9.  Also,  the  algorithm  correctly  estimates  the  pixel  bias  confirming  that  the  bias  must 
arise  from  a  noise  source  following  the  Poisson  distribution  (i.e.  dark  current). 


second  surface  pixel,  (32,18),  where  the  raw  waveform  results  in  an  incorrect  range  determina¬ 
tion.  In  contrast,  the  object  recovery  algorithm  (GEM0)  yields  an  improved  range  estimate  by 
sufficiently  estimating  the  true  waveform. 

5  CONCLUSIONS 

Utilizing  waveform  sampling  capability,  the  positive  effects  of  object  recovery  in  3D  FLASH 
LADAR  range  estimation  is  clearly  evident.  The  innovative  3D  FLASH  LADAR  sensor  pro¬ 
vides  both  an  imaging  and  ranging  ability  enabling  established  theory  to  be  applied  to  a  novel 
manner.  Given  simulation  and  experimental  results,  it  is  clear  the  chosen  model  and  noise 
sources  are  an  appropriate  choice  for  3D  FLASH  LADAR  data  operating  under  certain  con¬ 
ditions  (SULAR  mode  meeting  spatial  sampling  requirements).  The  raw  data  coming  off  the 
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sensor  does  not  fit  the  model,  but  straight-forward  pre-processing  steps  convert  the  data  to  an 
acceptable  form  for  the  algorithms. 

In  mild  spatial  blurring  conditions,  simulation  results  predict  that  the  GEM  algorithms  in¬ 
crease  range  estimation  performance  substantially  over  no-processing  and  the  Wiener  filter 
method.  Again,  the  Wiener  filter  even  has  an  unfair  advantage  because  it  is  provided  with 
the  exact  PSF  function  used  in  generating  the  data  while  the  GEM  algorithms  have  to  estimate 
the  PSF.  Considering  the  experimental  data,  its  performance  is  nearly  diffraction-limited  as  ev¬ 
idenced  by  the  experimental  PSF  and  OTF.  However,  the  GEM  algorithms  still  increase  range 
estimation  performance  over  the  Wiener  filter.  Supported  by  simulation  results,  it  is  appropriate 
to  say  that  the  GEM  algorithm  would  show  even  better  range  estimation  performance  versus  the 
Wiener  filter  in  severe  isoplanatic  atmospheric  blurring  conditions  or  with  sub-optimal  optics. 

A  trade-off  exists  for  Wiener  filter  and  object  recovery  algorithms  between  computation 
cost  and  range  accuracy.  The  Wiener  filter  is  the  least  computationally  taxing  object  recovery 
algorithm,  but  is  the  least  accurate  and  requires  a  priori  knowledge  of  the  PSF.  The  GEM  algo¬ 
rithms  are  computationally  expensive,  but  provides  the  best  range  performance  and  can  perform 
blind  deconvolution.  Considering  the  GEM  algorithms,  the  pulse-shape  estimator  is  extremely 
valuable  in  that  it  can  perform  range  estimation  on  single  cube  thereby  removing  potential  for 
any  registration  or  timing  errors.  If  multiple  cubes  are  available  and  properly  registered,  object 
estimation  is  undoubtedly  the  best  algorithm  to  use.  Although,  none  of  the  algorithms  were 
able  to  match  the  success  found  in  simulation.  Any  residual  error  in  the  experimental  results 
can  be  attributed  to  system  noise,  the  detected  light  containing  residual  laser  speckle,  residual 
gain  error,  and  detector  blurring. 

There  are  prospective  avenues  for  continued  investigation  and  improvement.  The  pulse- 
shape  estimation  is  very  dependent  on  the  selected  waveform  model.  Improvements  in  the 
range  estimation  would  be  realized  if  a  true  waveform  model  for  the  transmitted  laser  pulse  was 
derived  or  calculated  experimentally.  Errors  in  the  experimental  data  result  from  assuming  a 
generalized  shape  that  is  corrupted  by  distorting  effects  (spatial  blur,  pixel  blur,  and  noise).  In 
addition,  the  variable  of  interest  (range  term)  would  ideally  be  directly  estimated.  The  maximum 
likelihood  solution  for  the  range  term  could  be  achieved  if  another  model  was  discovered.  The 
algorithm  in  this  paper  extracts  the  range  from  the  maximum  likelihood  solution  for  the  pulse- 
shape.  Also,  even  after  the  pre-processing  steps,  the  experimental  data  exhibits  noisy  behavior. 
A  more  thorough  characterization  of  the  3D  FLASH  LADAR  noise  sources  would  augment  or 
verify  the  chosen  noise  sources.  Finally,  isoplanatic  imaging  is  valid  for  the  experimental  set¬ 
up  in  the  laboratory.  However,  object  recovery  from  3D  FLASH  LADAR  observations  subject 
to  heavy  anisoplanatic  turbulence  would  provide  an  ability  to  improve  range  estimation  in  a 
variety  of  field  or  operational  situations. 
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